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Abstract. Artificial viscosity is needed in Smooth Parti- 
cle Hydrodynamics to prevent interparticle penetration, to 
allow shocks to form and to damp post shock oscillations. 
Artificial viscosity may, however, lead to problems such as 
unwanted heating and unphysical solutions. A modifica- 
tion of the standard artificial viscosity recipe is proposed 
which reduces these problems. Some test cases discussed. 
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1. Introduction 

Smooth particle hydrodynamics (SPH) has become an in- 
creasingly popular method in simulations of different as- 
trophysical phenomena. Its Lagrangian formulation allows 
the study of large density differences. The particle formu- 
lation makes grids unnecessary and allows for easy im- 
plementation in three dimensions. A recent review of the 
method can be found in Monaghan (1992), and a technical 
description in for example Hernquist & Katz (1989). 
Two forms of artificial viscosity are needed in SPH, the 
bulk and the von Neumann-Richtmyer artificial viscosity 
respecticely. They prevent interparticle penetration, allow 
shocks to form and damp the post shock oscillations. They 
do, however, induce transfer of kinetic energy in fluid mo- 
tions to thermal energy. Many simulations using SPH in- 
volve the compression of a gas, often in a gravitational 
collapse of an initially cool gas cloud. The velocities in 
these simulations can be highly supersonic, which implies 
that the unphysically large artificial viscosity may domi- 
nate the heating. This heating and deceleration of the gas 
prevent further collapse. 

Martel et. al. (1995) have introduced a new formalism 
which they call Adaptive SPH. They use an ellipsoid ker- 
nel that adapts itself to the surrounding medium, and 



therefore avoid unnecessary use of artificial viscosity. The 
authors claim good results in cosmological collapse simu- 
lations. 

In the present study a different approach is suggested. The 
von Neumann-Richtmyer artificial viscosity is restricted 
to supersonic velocities. It is therefore used only to form 
shocks and to prevent interparticle penetration for super- 
sonic particles. Its adverse effect at subsonic velocities is 
avoided. A region under compression is described by a 
limited number of particles. The relative velocities among 
neighbouring particles in a region under compression can 
be supersonic due to the limited resolution, but despite 
that the gas is not expected to form shocks. Therefore the 
von Neumann-Richtmyer viscosity is restricted to regions 
that are not under compression. To avoid spurious heating 
in the subsonic regions, but prevent interparticle penetra- 
tion and maintain the damping of post shock oscillations, 
the bulk artificial viscosity is modified. Instead of integrat- 
ing the effect from the neighbouring particles individually, 
their collective effect at one point is considered. 
The SPH method is briefly described, and the suggested 
changes in artificial viscosity are presented. These modifi- 
cations of artificial viscosity have been tested in the sim- 
ulation of a shock tube, the homologous compression of 
a gas sphere and the gravitational collapse of an initially 
cool gas sphere at rest. The results are compared with 
results from a standard formulation of artificial viscosity. 



2. The SPH Methodology 

In SPH one models a number of particles that carry the 
physical quantities, where the particles' distribution in 
space describes the density distribution. To simulate a 
fluid each particle's mass is smoothed over a radius r. 
Following for example Hernquist & Katz (1989), such a 
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smoothed quantity at r be written for N particles as 



N 



/(r) = f dr'f{r)w{r' -r,h) = J2 ^firj)w{rij,hij)il) 



where niij + (m,; + m.j)/2, p^j = {p, + Pi)/2. h^j = (/?.,; + 
hj)/2 and — Vi — Vj. The smoothing kernel, w, has the 
property 



j> drw{r, h) 



1. 



(2) 



The kernel used here is the spline kernel from Monaghan 
and Lattanzio (1985): 



wir 




1<X<2 



^ h 

h — 



(3) 



This is a smooth kernel with compact support over the 
radius 2h around the particle. The smoothing length h 
is varying in time and updated every iteration to keep 
the interactions with other particles to a specific number. 
They are called the particle's neighbours, and the number 
of neighbours for each particle in the tests in this paper is 
64. 

In the SPH formulation there are different forms of the 
discretization of the Navier Stokes equations. Here the ex- 
pressions of Hernquist & Katz (1989) are used, i.e. 



dVi_ ^ 
dt 

duj 

dt ~ 



2 1^3 



rrii 



(4) 



where Vjj = — Vj and = {Pi + Pj)/2. The continuity 

equation is automatically satisfied due to the Lagrangian 
formulation, and the density calculated from Eq. (1) be- 
comes 



Pi 



N 



(5) 



The standard artificial viscosity term, Ily , is defined 



P'ij 



_ , 1 if r,j 

"'^"10 if Tij 
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(6) 



where = (cj + Cj)/2. The first and second term in the 
expression of Iljj in Eq. (6) represents the bulk and the 
von Neumann Richtmyer artificial viscosity respectively. 
The constant v = 0.01ft, is a fudge parameter to prevent 
the artificial viscosity to become too large. The artificial 
viscosity is used only when Tij ■ Vij < 0, that is when two 



particles are approaching each other. To close the system, 
the pressure is defined as P + (7 — l)pu, where u is the 
thermal energy density and 7 the adiabatic index. The 
particles' quantities are updated using a standard leapfrog 
integrator with the time step 



dt = Mm[5ti 



ehi 



Vi + Ci + l.2{aCi + (iHi^max) ' 



(7) 



where e = 0.3 is the Courant factor to stabilize the inte- 
gration, and Umax the maximum n from the interactions 
with the other particles. Since all particles are integrated 
with the same time step, the smallest time step from all 
particles is used in the integration. In the leapfrog integra- 
tor the velocity and internal energy density are integrated 



at half time steps, t„_i/2,i 



-1/2; 



while the position is 



integrated at whole time steps, tn,tn+i, as 



1-1-1/2 



= V,- 



-1/2 ■ 



St 



dWi 
dt ,„ 
5t 



^i,ra-|-l/2 — Wi,n-l/2 i TJj , 



(8) 




The viscous acceleration terms in Eq. (4) scales with a = 
/3 = 1 as 



(9) 



The artificial viscosity can therefore lead to undesirable 
effects, because the velocity differences are smoothed on 
a time scale of roughly h/{c + v). The velocities in the 
model will therefore be smoothed out unless it expands or 
if there is some driving mechanism such as gravitation. To 
conserve the energy the particles are heated, which may be 
unphysical. The heated gas may reach an equilibrium state 
earlier than expected. A way of preventing interparticle 
penetration without unnecessary heating of the gas could 
therefore be useful, and this implies that there is a need 
to restrict the artificial viscosity to the shocks as much as 
possible. 

3. Restrictive Use of von Neumann-Richtmyer Ar- 
tificial Viscosity 

In SPH the formation of shocks is mostly an effect of the 
von Neumann-Richtmeyer viscosity. To avoid the unde- 
sirable effects in the subsonic region, I propose to restict 
the use of it to those regions. This modification will allow 
shocks to form, and prevent interparticle penetration at 
supersonic velocities. 

Consider a gas cloud in a spherical symmetric gravita- 
tional collapse. Suppose the physical viscosity is small, 
so that the compression can be assumed to be adiabatic. 
When the model has reached an equilibrium, the pres- 
sure force that prevents further gravitational compression 
balance the gravitational force. If the cloud is warm, the 
forces may balance each other even in its initial state. But 
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if on the other hand the cloud is initially cool, it must 
be compressed to gain the required thermal energy den- 
sity, perhaps by orders of magnitude. In many cases this 
is not possible with standard artificial viscosity, Eq. (6), 
due to poor resolution. If the cloud in the example above 
has a radius of unit length and is modelled with 10'^ par- 
ticles, the mean interparticle distances are about 0.1. In 
the spherical compression the particles at different radii 
therefore have supersonic relative velocities if the gas is 
cool enough. Standard artificial viscosity, Eq. (6), will de- 
celerate and heat the gas, and prevent further compres- 
sion. 

Since p oc the time derivative for any particle can be 
written 



P = — 



3ph 



(10) 



This relation can be used to decide whether a particle fol- 
lows the fluid or if artificial viscosity is necessary. Consider 
two particles with a separation of r moving towards each 
other with a speed of r they follow the fluid in the neigh- 
bourhood, the neighbourhood is under compression and 
the particles' relative velocity will approximately satisfy 



h 
h' 



(11) 



If two particles are approaching each other at a velocity 
exceeding the sound speed, that is if 



Tj,- • Wij < and va > Ci- 



(12) 



artificial viscosity is necessary to prevent interparticle pen- 
etration. I therefore propose to restrict the von Neumann- 
Richtmyer artificial viscosity as 



i-^vNR.ij — -^r— 

■ ■' Ptj 

Mij — iT-ij r'-'.+n'-f. "ij 
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(13) 



4. Modified Bulk Viscosity 

A smoothed quantity can be calculated at any point in 
the fluid by Eq. (1). The smoothed velocity at r, is 



Af N _ 



Pij 



(14) 



If this point coincides with a position for a particle, the 
smoothed and the particle's individual velocity will be 
different in general. This difference is used to construct 
a modified bulk viscosity. Benz (1990) suggests that this 
quantity could be used when integrating the position to 



prevent interparticle penetration. It is true that it prevents 

interparticle penetration, but it unfortunately introdiices 
conservation problems when the particles are not moved 
at their individal velocity. If the contribution from parti- 
cle i subtracted, it does however say something about the 
fiuid around the particle. The smoothed velocity at r can 
then be redefined as 



Ar' = 




(15) 



where Wij = w(rij/hij = 0) = l/nh'^j. The denominator 
scales the expression so that Avj = Vj, if for all j vj = Vj. 
In the standard formulation of the bulk artificial viscosity 
particle i interacts with each neighbour separately, where 
the acceleration and time derivative of the internal en- 
ergy is added to the particle as described in Sect. 2. This 
introduces problems described in Sect. 3. The individual 
velocity of particle i can then be seen as a deviation from 
the smoothed value at the particle's position and the arti- 
ficial viscosity as a correction to the individual velocity. I 
propose a modified bulk viscosity to replace the standard 
bulk viscosity, where the fiuid around the particle is con- 
sidered from a collective contribution from the neighbours 
in one single interaction. This modified bulk artificial vis- 
cosity is defined as 



Si 



dt ~ hi 

1 if V Avi < 
if V • Avi > 



(16) 



where Nn the number of neighbours. The constant r] 

around unity, and used in the same way as the constants a 
and /? in Eq. (6). This interaction can be seen as if the par- 
ticle interacts with a virtual particle with the smoothed 
velocity Avj, have a mass of N„m,i and lies at a distance of 
h the direction of Vj . This expression thus becomes similar 
to Eq.(6) and affects the same velocity regime. The dif- 
ference is that in Eq. (16) the collective contribution from 
all neighbours is considered in one single interaction. 
Now consider the integration of the velocitites and inter- 
nal energy density from time i„_i/2 to i„+i/2 with the 
time step 5t. From Eq. (8) the velocity for particle i at 
tn+i/2 is 



Vi,n-|-l/2 — Vi_„_i/2 + 

which gives the change in kinetic energy for particle i: 



(17) 



AI^, 



1-1-1/2 



dt 



2 dt 



-1/2 



dt 



■Si 



r~-(i8) 



If AUi = rrii {dui/dt)St is the change in internal energy 
for the particle it is possible to conserve the energy, that 
is require that AT, -|- MJi = 0. The time derivative of 
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the internal energy for the particle is then consequently 
defined as 



dt 



-Vi n-1/2 ; ; ; Ot (19) 

' dt 2 dt dt ^ ' 



to conserve the total energy. 

The artificial also must prevent particle penetration. The 
modified artificial viscosity is calculated with respect to 
an integrated mean of the neighbours. Therefore a particle 
will move less than eh from the definition of the time step, 
Eq. (7), regardless of the neighbours' individual sound 
speed and individual velocities. This should be compared 
with the distance to the closest which are approximately 
one h with 64 neighbours. Since the neighbours have no 
identity, this may lead to penetration with a few parti- 
cles which have a sufiicient deviation from the integrated 
mean. 

If the particles i and j are each others neighbours and 
that their other neighbours give the same contribution to 
their respective velocities, one concludes from Eq. (15) and 
(16) that the impulse is conserved. Their other neighbours 
do, however, not give the same contribution, duo to the 
limited resolution. Tests of self gravitating rotating disks 
show that the angular momentum and impulse are well 
conserved. 



^/-directions. At the ends of the tube, z = —5 and z = 4, 

no boundary conditions were appUcd, so the particles are 
allowed to move away from the tube. Their velocities are 
however low compared with the velocity of the shock and 
do not affect the shock model. The 16384 particles are 
ordered in a cubic centered grid, such that 8192 particles 
are distributed at— 5<^;<— 4, and the rest at — 4 < 
z < 4. With a total mass of 1.0 mass units the density 
distribution is 



_ ri/4 where 
^ \l/16 where 



-5 < < -4 
-4 < 2; < 4 



(20) 



The initial thermal energy density is set to 0.01, and the 
particles have no initial velocity. A shock is formed at the 
discontinuity at 2; = —4, which propagates to the right. 
This is a rather weak shock, which is a better test than a 
strong shock. The reason is that here the particle veloc- 
ities are not much larger than the sound speed, because 
if they are the standard and modifies artificial viscosity 
become similar. Fig. 1 shows the shocks at f = 20 with 
the standard artificial viscosity, Eq. (6), and Fig. 2 using 
modified artificial viscosity, Eq. (13) and (16). A compari- 
son between Fig. 1 and 2 shows that the shocks formed by 
the two versions of artificial viscosity are similar, so that 
the modified artificial viscosity is able to work in the same 
way as the standard artificial viscosity. 



5. Tests 

Any form of artificial viscosity must be able to form and 
propagate a shock. The modification of the artificial vis- 
cosity, Eq. (13) and (16), is tested in a shock forming test 
and compared with the standard artificial viscosity, Eq. 
(6), introduced by Monaghan and Gingold (1983). The 
ability of the modified artificial viscosity to compress the 
gas without viscous deceleration and heating has also been 
tested in a homologous compression of a gas sphere. The 
test constructed by Evrard (1988) is used to study the dif- 
ferences between the artificial viscosities in a gravitational 
collapse of an initially cool gas cloud. 
The number of particles is varied between 8192 and 16384, 
and the number of neighbours for each particle is 64, so 
that h is varying in time and space. In the equation of 
state the adiabatic index is 7 = 5/3 to model an ideal 
monoatomic gas. Dimensionless units are uscid to keep the 
quantities in the model around unity, where the gravita- 
tional constant G = 1. In a model the total mass M = 1, 
the typical length L = 1 and time T — I, which relates 
to the gravitational constant as G = L^/MT"^. The real 
quantities of the model can be calculated by inserting the 
corresponding quantities in the desired unit system. 



5.1. Shock Formation Test 

A box with dimensions dx x dy x dz = 1 x 1 x 9 is used. 
Periodic boundary conditions are applied in the x- and 



Fig. 3. The energy curves for the homologous compression of 
a gas with different forms of artificial viscosity. The solid line 
represents the compression without any artificial viscosity at 
all, the dashed the modified artificial viscosity, Eq. (13) and 
(16), and the dotted the standard artificial viscosity, Eq. (6). 
The curves that decline in the beginning represent the kinetic 
energies for the three models, those that rise represent the 
internal and the uppermost straight curves represent the total 
energy. 
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Fig. 1. A shock in a shock tube with standard artificial viscosity, Eq. (6), which was formed by the density discontinuity 
described in Eq. (20). Fig. la shows at time t = 20 the velocity distribution, lb the pressure and Ic the density distribution. 



Fig. 2. Same as Fig. 1, but with modified artificial viscosity, Eq. (13) and (16). 



Fig. 4. The density distribution for the homologous compression of the sphere at maximum compression without any artificial 
viscosity in 4a, with standard artificial viscosity in b, Eq. (6), and modified artificial viscosity, Eq. (13) and (16) in 4c. The 
results in Fig. 4a and 4c are plotted with the same scale, while another scale must be used in 4b. The solid line in Fig. 4a and 
c represent the zeroth order theoretical estimate from Eq. (26) at maximum compression where the radius is 0.029 and the 
density 10000. 



5.2. Homologous Compression of a Gas Sphere 

The ability of the modified artificial viscosity, Eq. (13) and 
(16), to handle compression of a gas realistically is tested. 
Initially 8192 particles are distributed unifomly in a sphere 
with radius Rinit = 1 on a shghtly disturbed cubic cen- 
tered grid. There are no boundary conditions applied, but 
there is an initial velocity distribution directed towards 
the origin according to 

^init{r) = -Vo-^, (21) 



where Vq = 2 is a constant and r the position in the sphere. 
The particles arc initially isothermal with a thermal en- 
ergy density of Umit = 0.001, and with a total mass of 
M = 1 , which gives an initial density of Pmit = 3/47r. The 
sum of the total kinetic, Ekin, and thermal, Eth, energies 
is 

Etot = Ekin + Eth = ^dm + uM = ,^2) 
/(f ^inr'^pdr + /xM = | + 0.001 = 1.201 ^ ^ 

Assume that the compression is adiabatic, so that Pois- 
son's equation, 

P = Kp'^, (23) 
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is valid. Prom the equation of state, P = (7 • 
the adiabatic index, 7 = 5/3, it follows that 



1.73-10- 



l)up, and 



(24) 



If it is assumed that all kinetic energy is converted to ther- 
mal energy at maximum compression, the thermal energy 
density at this point is 



M 



1.201. 



(25) 



This gives a density and radius of the compressed gas 
sphere as 



P = 



K 



(7-1)1. 



R' 



( 3M A ^ _ 
\inp' ) - 



(26) 



0.029 



This zeroth order approximation is useful to compare with 

calciilations with different forms of artificial viscosity. The 
initial conditions also have the advantage that no artificial 
viscosity is needed to prevent interparticle penetration. A 
small initial pressure is sufficient to dcccllcratc the parti- 
cles to zero and prevent them to move through the origin. 
The results from the test with modified artificial viscosity 
can therefore not only be compared with the analytical ap- 
proximation Eq. (25), but also with a model without any 
artificial viscosity. The energy curves from such a compar- 
ison are shown in Fig. 3. In Fig. 4 the density distributions 
at time t = 0.5 are compared. 

The heating in the case with standard artificial viscosity, 

Eq. (6), starts immediately, because it depends on the rel- 
ative velocities, while the heating without artificial viscos- 
ity and with the modified artificial viscosity is negligable 
until t « 0.3. The modified artificial viscosity has a little 
less steep energy curve compared with the case without 
artificial viscosity, but gives an almost a compressed gas 
as without artificial viscosity as seen in Fig. 4. The true 
density distribution is not known, but the zeroth order ap- 
proximation from Eq. (25) gives approximately the same 
size as the model without artificial viscosity. 



5.3. The Evrard Gravitational Collapse 

The initial conditions arc those from Evrard (1988), which 
is a isotherm sphere at rest with a thermal energy density 
of u = 0.05, radius R = 1 and mass M = 1. Initially 8192 
particles arc distributed uniformly in a sphere with radius 
i? = 1 on a slightly disturbed cubic centered grid with a 
density distribution of 



p{r) 



M 1 



2TTRr 



(27) 



Standard, Eq. (6), and modified artificial viscosity, Eq. 
(13) and (16) are tested with this model and compared 



Fig. 5. The energy curves for the Evrard gravitational collapse 
using standard artificial viscosity, Eq. (6), represented by a 
dashed line, and modified artificial viscosity, Eq. (13) and (16), 
represented by a solid line. The uppermost curves represent the 
internal thermal energy, the next the kinetic, the straight line 
the total and the two curves below the others' the potential 
energy. The solid line is from an accurate one-dimensional PPM 
calculation at t = 0.77 from Steinmetz & Miiller (1993). 



with each other. 

Due to the cold initial state the sphere begins to collapse. 
Since the central part is more dense than the outer payers, 
the collapse is more rapid around the origin. A high cen- 
tral pressure and density is build up, and the central parts 
starts to expand at t ~ 0.8. Where the expanding parts 
meet the infalling gas, an outward propagating shock front 
forms. Eventually the gas reach a virial equilibrium. 
The total energies are shown in Fig. 5a and b. The curves 
are more shallow in Fig. 5a, where standard viscosity was 
used, than in Fig. 5b. In Fig. 6 the velocity, density and 
pressure distributions are plotted with standard artificial 
viscosity at t = 0.8 when the shock is formed. This can 
be compared with the results from modified artificial vis- 
cosity for the same quantities which are shown in Fig. 7. 
The main difference between these models are the sharper 
gradients with modified artificial viscosity. This is an ef- 
fect of the ability to compress the gas with modified arti- 
fial viscosity. The infalling gas is allowed to move inwards 
without decelleration until it meets the shock front. 



6. Conclusions 

In SPH two forms of artificial viscosity is necessary, the 
bulk and von Neumann-Richtmyer viscosity. The standard 
recipy does however induce smoothening of the velocity 
differences and heating of the gas. Another problem is 
that standard artificial viscosity prevents compression of 
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Fig. 6. The velocity, density and pressure distributions for the Evrard gravitational collapse using standard artificial viscosity, 
Eq. (6). The radial velocity, density and pressure arc plotted at t=0.8. The crosses are values at t=0.77 are results from an 
accurate one-dimensional PPM calculation from Steinmetz & Miiller (1992). 



Fig. 7. Same as Fig. 6, but with modified artificial viscosity, Eq. (13) and (16). 



a gas. 

To partly overcome some of these problems, I propose a 
modification of the artificial viscosity. The von Neumann- 
Richtmycr artificial viscosity is used only in supersonic 
regions where the gas is not under compression. The bulk 
viscosity is modified to be calculated from an integrated 
mean of the velocities in the neighbourhood of the parti- 
cle, instead of calculating the viscous interaction from each 
neighbour separately. Tests cases have been performed to 
test these abilititcs of the proposed modified form of artifi- 
cial viscosity, and to compare them with standard artificial 
viscosity. 
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